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Abstract. The shear strength of cohesionless granular materials is generally attributed to the com- 
pactness or anisotropy of their microstructure. An open issue is how such compact or anisotropic 
microstructures, and thus the shear strength, depend on the particle properties. We first recall the 
role of fabric and force anisotropics with respect to the critical-state shear stress. Then, a model of 
accessible geometrical states in terms of particle connectivity and contact anisotropy is presented. 
This model incorporates in a simple way the fact that, due to steric exclusions, the highest levels 
^ of connectivity and anisotropy cannot be reached simultaneously, a property that affects seriously 

i the shear strength. We also analyze the force anisotropy in the light of the specific role of weak 

forces in sustaining strong force chains and thus the main mechanism that underlies anisotropic 
force patterns. Finally, we briefly discuss the effect of interparticle friction, particle shape, size 
polydispersity and adhesion. 
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3 ■ 1 Introduction 

oo ' 

f^) , Since the early work of Coulomb in 1773, the plastic yield behavior of granular materials 

has remained an active research field in close connection with soil mechanics and pow- 
der technology [Mitchell and Soga(2005)| |Nedderman(l 992)1 . According to the Mohr- 



X 



Coulomb yield criterion, for normal and shear stresses a and r acting on a slip plane, the 
plastic threshold r c is the sum of two terms: 

r c = c + a tan tp, (1) 

where c is a cohesive strength and <p is the internal angle of friction depending only on the 
nature of the granular material. This criterion expresses the pressure dependence of shear 
strength which is a distinctive feature of granular media. Given ([T), the shear strength 
of cohesionless materials (c = 0) can be represented by the (dimensionless) stress ratio 
T c /a = fj, c = tan p. Since the angle ip is a bulk property, it can be expressed in terms 
of stress invariants. Let a a (a = 1, 2, 3) be stress principal values. The average stress 
is p — (<ti + <72)/2 in 2D and p — (ai + <72 + 03)/3 in 3D. We define the stress 



deviator by q = (cri — <72)/2 in 2D and q = {a\ — cr^)/Z in 3D under axisymmetric 
conditions (02 = 0-3). With these notations, it can be shown that sin ip = q/p in 2D and 
sin if = 3q/ (2p + q) in 3D. 

This picture of shear strength in granular media holds as a basic fact although the complex 
plastic behavior of granular media can not be reduced to a single strength parameter. In 
particular, the shear strength and plastic flow (dilatancy) depend on the granular structure 
and direction of loading, the latter reflecting the anisotropy of the structure. Since the 
shear strength is state-dependent, it cannot be considered as a material property unless 
attributed to a well-defined granular state. The internal angle of friction ip is often associ- 
ated with the critical state (steady state or residual state) reached after long monotonous 
shearing; see Fig. Q] This state is characterized by a solid fraction p c independent of the 
loading history and initial conditions |Wood(1990)| . 

The critical-state strength is below the peak shear stresses occurring for dense states with 
solid fraction po > p c , but these states are metastable and often lead to strain localization 
|Darve and Lao uafa(2000)| |Vardoulakis and Sulem(19 95)|. For loose states with po < 
p c , the critical state is reached asymptotically following diffuse rearrangements. Hence, 
apart from these transients, which are governed by the evolution of internal variables 
pertaining to the microstructure and are important in formulating elasto-plastic models, 
the critical-state shear strength represents a stable plastic threshold for a granular material. 
In this paper, we are interested in the critical-state strength as a material property of 
cohesionless granular materials. The critical-state friction angle ip c can be described as a 
coarse-grained (or homogenized) friction angle between two granular layers sliding past 
each other. Nevertheless, the macroscopic status of ip c as a Coulomb friction angle, on 
the same grounds as those of dry friction between solid bodies, should not eclipse the 




Figure 1 : Normalized shear stress as a function of cumulative shear strain in a 2D simple 
shear simulation by the contact dynamics method for two different values of the initial 
solid fraction. 



fact that the granular friction angle is a bulk property to which adequate tensorial stress 
analysis should be applied (this was indeed the contribution of Mohr) and where the slip 
planes are not a priori defined, in contrast to solid friction which is a surface property 
at the macroscopic scale [Radjai e t al.(2004))|. Depending on the boundary conditions, 
the critical state occurs either homogeneously in the whole volume of a granular sample 
or inside a thick layer of several particle diameters in the advent of strain localization 
IBardet and Proubet( 1992)1 |Herrmann et al.(1995)| |Vermeer( 1990)1 |Moreau( 1997| . In 
both configurations, tp c stems from various granular phenomena such as friction between 
particles, anisotropy of the microstructure, organization of force networks and dissipation 
due to inelastic collisions. We consider below these effects and their respective roles in 
enhancing or restraining granular friction. 



2 Effect of interparticle friction 

While solid friction between particles underlies the frictional behavior of granular mate- 
rials, it is not obvious how and through which physical mechanisms it comes into play. If 
shear deformation took place as a result of sliding between all contacts along a slip plane, 
the friction angle ip c would simply echo the friction between particles. An example of 
such a configuration is a regular pile of cubic blocs subjected to a vertical load. Horizon- 
tal shearing of this pile implies sliding between at least two rows so that the shear strength 
of the pile is a straightforward effect of friction between the blocs. However, discrete nu- 
merical simulations suggest that in sheared granular materials, rolling prevails over sliding 
[Radjai era l.(1998)|. In quasistatic shear, sliding occurs at only ~ 10% of contacts, and 
these sliding contacts belong essentially to weak contacts (see below) oriented on average 



along the minor principal stress directions | Radjai et al.(1999) Staron and Radjai(2005) 



|Staron et al.(2005)) . Hence, the relationship between (p c and the local friction angle ip s 
involves the inhomogeneous distribution of forces and mobilization (or activation) of the 
friction force at rolling contacts. 

This relationship is far from linear as shown in Fig. |2] The critical-state coefficient 
/i c = tan ip c is above /i s ~ tan tp s at small values of the latter, and at larger values it tends 
to a plateau /i^ < /j, s [Corriveau et al.(1997)[|Tab oada et al.(2006)| . The transition from 
/i c — ii s < to yu c — fi s > occurs at /x c = p s ~ 0.5. Beyond fi s = 0.5, fi c is practically 
independent of [i s . The independence of <p c with respect to ip s at large values of the latter 
indicates that the role of interparticle friction is more subtle than expected from simple 
models. Moreover, the nonzero value of ipo shows clearly that the interparticle friction is 
not the only source of frictional behavior in the critical state | Roux and Radjai(2001)) . 



The direct contribution of interparticle friction to shear strength, i.e. without interposi- 
tion by the microstructure as will be analyzed below, may be evaluated from a decom- 
position of the shear stress. The stress tensor a a p in a control volume V can be ex- 
pressed as |Rothenburg and Selvadurai(1981)||Christoffersen et al.(1981)||Moreau(1997)| 



Bagi(1999)||Staron et al.(2005jr 

<J*p = n t <4#>, (2) 



where rif, is the number density of bonds (contacts), t a is the a-component of the branch 
vector t joining the centers of particles at contact i and fl is the /3-component of the 
force vector / acting at the contact i between the two particles. 

The contribution of friction forces can be estimated by replacing in equation © / by 
f-tt, where t is the unit vector along the friction force. The contribution of normal forces 
is the complementary tensor obtained by replacing / by / • n n, where n is the unit vector 
perpendicular to the contact plane. The corresponding shear strengths q t and q n can then 
be calculated in the critical state. Numerical simulations show that the ratio q t /q is quite 
low (below 10%) |Cambou( 1 993)[ . This counterintuitive finding underlines the role of 
interparticle friction as a parameter acting "behind the scenes" rather than a direct actor 
of shear strength. Our simulations show that, due to disorder and force/moment balance 
conditions as well as kinematic constraints such as rotation frustration, the friction forces 
inside a granular packing are strongly coupled with normal forces. For example, highly 
mobilized friction forces are rare events and the distribution of friction forces reflects for 
the most part that of normal forces. We consider below such effects in connection with 
granular microstructure. 

3 Harmonic representation of the microstructure 

The microscopic expression of the stress tensor in equation (0 is an arithmetic mean in- 
volving the branch vectors and contact forces. Hence, for analyzing the particle-scale 
origins of the shear strength, we need a statistical description of the granular microstruc- 
ture and force transmission. Noticing that the shear stress corresponds to the deviation 
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Figure 2: The critical-state friction coefficient fj, c as a function of sliding friction coeffi- 
cient /i s between particles in biaxial shearing of a sample of 5000 particles. 




of stress components from the mean stress p = tr(cr)/d (for space dimension d) along 
different space directions, the useful information for this analysis is the density and av- 
erage force of all contacts pointing in the same direction as a function of this direction. 
These functions can be expanded in Fourier series in 2D and in spherical harmonics in 



3D[Rothenburg and Bathurst(1989) Ouadfel and Rofhenburg(2001)|. Since the contacts 



have no polarity, the period is n. 

For illustration, we consider here only the 2D expansions truncated beyond the second 
term: 

P e {6) = ±{l + acos2{6-9 b )}, 

(f n ){0) = (/){l + a„cos2(0-0 n )}, (3) 
(ft)(0) = </)a t sin2(0-0 t ), 

where Po is the probability density function of contact normals, and /„ and ft are the force 
components along (radial) and perpendicular to (orthoradial) the branch vector, respec- 
tively. The parameters a, a n and at are the anisotropies of branch vectors, radial forces 
and orthoradial forces, respectively, 9 b , 6 n and 8 t being the corresponding privileged di- 
rections. The sine function for the expansion of the orthoradial component f t is imposed 
by the requirement that the mean orthoradial force is zero to satisfy the balance of force 
moments over particles whereas the mean radial force {/} is positive (repulsive). We also 
note that for circular and spherical particles the radial and orthoradial force components 
coincide with normal and tangential forces, respectively. 

This harmonic representation with only three anisotropy parameters provides a good ap- 
proximation for numerical data. Using the functions (0), the stress components a a j3 can 
be written as an integral over space directions: 

<J a0 = n b (£) {</„)(0)n„(0) + (f t )(9)t p (9)} P e (9) dS, (4) 
Jo 

where n x = cos(8) and n y = sin(0), t x = — sin(0) and t y = cos(#). It has been also 
assumed that the branch vector lengths £ are not correlated with forces. 
Equation (j4]i together with the harmonic approximation expressed in equation ([3} yield 
the following expression for the normalized stress deviator | Radjai et al.(2004)) : 

^ ~ I {a cos 2{0 a - 6 b ) + a n cos 2(0 a - 6 n ) + a t cos 2{0 a - t )} , (5) 
p 2 

where 9 a is the major principal direction of the stress tensor. In deriving equation ©, the 
cross products among the anisotropies have been neglected. In the critical state, the privi- 



leged directions coincide, i.e. 9 b ~ 8„ ~ 9 t ~ 9 a , so that [Rothenburg and Bathurst(1989) 
Ouadfel and Rothenburg(2001)) 

— - o + a nc + a tc} , (6) 

P 2 

where the anisotropy parameters refer to the critical state. In 3D, a similar relation can be 
established by means of spherical harmonics (A zema et al. (2008)1 : 

— — 7 {a c + a nc + a tc ] (7) 
p 5 



These relations exhibit two microscopic sources of the shear strength in a granular pack- 
ing: 1) fabric anisotropy, represented by the parameter a and 2) force anisotropy, captured 
into the parameters a n and a*. Hence, the material parameters influence the shear strength 
via fabric and force anisotropics. For example, the saturation of <p c for <p 3 > 0.5 (sec- 
tion|2]means that, increasing the interparticle friction beyond this limit does not enhance 
anisotropy. 



4 Accessible geometrical states 

In this section, we focus on the fabric anisotropy a which represents the excess and loss 
of contacts along different space directions with respect to the average contact density. 
The latter is commonly represented by the coordination number z (mean number of con- 
tacts per particle). In a granular material, z is bounded between two limits z m i n and 
z max . The lower bound z min is dictated by the force balance requirement. For example, 
stable particles often involve more than three contacts in 2D and more than four con- 
tacts in 3D. On the other hand, the upper bound z max is constrained by steric exclusions 
|Troadec et al.(2002)]. For example, in 2D for a system of monodisperse particles, a par- 
ticle can not have more than 6 contacts. In practice, this limit is reduced to 4 as a result of 
disorder. 

Within the harmonic approximation, the geometrical state of a granular system is defined 
by its position in the space of coordinates (z, a). We define two limit states: 1) the loosest 
isotropic state, characterized by [z = z m i n ,a = 0), and 2) the densest isotropic state, 
characterized by (z = z max ,a = 0). These states can be reached only by complex 
loading. For example, it is generally difficult to bring a granular system towards a dense 
isotropic state via isotropic compaction. The reason is that the rearrangements occur 
mainly in the presence of shearing, and the latter leads to fabric anisotropy. 
It is natural to assume that all accessible geometrical states are enclosed between the two 
isotropic limit states. In order to represent the geometrical states, it is useful to define the 
state function 

E{6) = zP g {6) = ~{1 + acos2(6» - 9 b )}. (8) 

7T 

The two limit isotropic states are E min = z min /n and E max — z max /n. The assumption 
that the geometrical states are constrained to stay between the two isotropic limit states, 
implies that the anisotropy a can not exceed a maximum a max depending on the value of 
z. With harmonic approximation (8), we obtain 

a max (z) = min {2 (l - *f) , 2 {^f- l) } . (9) 

This function is shown in Fig. [3] By construction, a max (z mm ) = a max (z max ) = 0. The 
largest anisotropy is 

a c = a m ax{z c ) = 2 — , (10) 

^max 1 drain 

with z c = (z m i n + z max )/2. According to equation (TTOb , a max increases with z for 
z < z c , and it declines with z for z < z c . When a = a c is reached along a monotonic path, 



neither anisotropy nor coordination number evolve since both contact gain and contact 
loss are saturated. In this picture, the critical state corresponds to the intersection between 
the two regimes with z = z c and a = a c . In 2D with weakly poly disperse circular 
particles and /u s > 0.5, a good fit is provided by assuming z m in = 3 and z rnax = 4. This 
yields z c = 3.5 and a c = 2/7. For lower values of fj, s , a c declines. 




Figure 3: Domain of accessible geometrical states based on the harmonic representation 
of granular microstructure. 
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Figure 4: Evolution of the geometrical state of a sheared packing for two different initial 
states simulated by the contact dynamics method. 



Fig. H] shows the evolution of a with z in simulated biaxial compression of two initially 
isotropic samples with initial coordination numbers zq = 3.1 and zq = 3.7. In both simu- 
lations, z tends to the same critical-state value z c ~ 3.35 with a c ~ 0.24. Remarkably, the 
anisotropy of the dense packing reaches and then follows closely the limit states. Equa- 
tion (O provides here an excellent fit to the data with only one fitting parameter z max . In 
the loose case, the trajectory remains entirely inside the domain of accessible states and 
the limit states are reached only at the critical state 

Equation ( TTOb predicts that the critical state anisotropy a c increases with z max — z m i n . 
The shape, size and frictional characteristics of the particles may therefore influence a c via 
z m i n and z max . For example, increasing the sliding friction between the particles allows 
for lower values of z min (stable configurations with less contacts) without changing z max 
(which depends only on steric exclusions) and leads to larger values of a c . 
One interesting aspect of the model of accessible states presented in this section is to show 
that the largest values of a and z can not be reached simultaneously. The critical value 
a c is not obtained with z max but with z c which is below z max . But higher levels of force 
anisotropics a nc and ate can be achieved with higher values of z. 

5 Weak and strong force networks 

According to equation ( fTOb , the shear strength is proportional to force anisotropics a nc 
and ate in the critical state. As for fabric anisotropy a, which was discussed in the last 
section, we would like to analyze here the mechanisms that underly force anisotropics. A 
basic feature of force distribution in granular media is the occurrence of numerous weak 
forces together with a subnetwork of strong forces appearing often sequentially (force 
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Figure 5: The partial shear stress, normalized by the mean stress, as a function of force 
threshold £. 



chains). The probability density function (pdf) P n (f n ) of normal forces in a macro- 
scopically homogeneous system in the critical state is such that more than 58% of con- 
tact forces are below the mean force (/„) and they have a nearly uniform distribution 
IRadjai et al.(1996)|[Mueth et al.(1998)l|Tsoungui et al.(1998)[ . These weak forces con- 
tribute only ~ 29% to the mean stress p. The pdf of strong forces (above the mean normal 
force (/„)) decays exponentially Radjail996a, Coppersmith 1996a, Radjail999, Majmu- 
dar2005, Metzger2004. The very large number of weak forces, reflecting the arching 
effect, is a source of weakness for the system. Weak regions inside a packing correspond 
to locally weak pressures and they are more susceptible to fail. A quantitative analysis of 
grain rearrangements indicates that during a quasistatic evolution those weak regions un- 
dergo local rearrangements, and nearly all sliding contacts are localized in weak regions 
IStaron et al.(2002)|[Staron et al.(2005)l|Nicot and Darve(2006)| . 

Let <S(£) be the set of contacts with a normal force /„ < £(/ n )- The set 5(oo) is the 
whole contact set. The weak and strong sets are <S(1) and <S(oo) — S(l), respectively. 
The partial shear stress q(£)/p and the fabric and force anisotropics a(£), a„(£) and at(£) 
can be calculated as a function of £ ]Radjai et al.(19 98)|. Our simulations show that 



q(£) ~ 0; see Fig. [5] This means that nearly the whole stress deviator is carried by the 
strong contact network, the weak contacts contributing only to the mean stress. Hence, 
the total stress tensor u is a sum of two terms: 

<T =p w I + <T S , (11) 




Figure 6: Weak and strong normal forces represented in two different grey levels. Line 
thickness is proportional to the normal force. 



where J is the unit tensor, p w is the weak pressure, and er s represents the strong stress ten- 
sor. Hence, from the stress transmission viewpoint, the weak contact set is a "liquidlike" 
phase whereas the strong contact set appears as a "solidlike" backbone transmitting shear 
stresses. The weak and strong networks are shown in Fig. |6]in thickness of segments join- 
ing particle centers for an assembly of 4000 particles subjected to biaxial compression. 
The zero shear stress in the weak network implies that, according to equation ©, at least 
one of the corresponding anisotropics is negative. Since the critical-state angles are as- 
sumed to be equal (6^ ~ 6 n ~ 9 t ~ 9 a ), a negative value corresponds to a rotation n/2 of 
the principal axes. Indeed, our numerical data show that the privileged direction of weak 
contacts is perpendicular to the major principal stress direction [Radjai et al.(1998)|. The 
strong forces occur at contacts that are, on average, aligned with the major principal di- 
rection of the stress tensor. Lateral weak forces prop the particles against deviations from 
alignment at strong contacts. In other words, the weak contacts play the same stabilizing 
role with respect to the particles sustaining strong forces as the counterforts with respect 
to an architectural arch. This bimodal transmission of shear stresses corresponds thus to 
a statistical description of arching effect in granular media. 

This stress-fabric correlation can be interpreted as a way for a granular system to optimize 
the shear strength. Indeed, the stress deviator q increases if a larger number of strong 
forces occur at contacts aligned with the major principal direction, implying thus a surplus 
of weak contacts in the perpendicular direction. This weakening of forces at contacts 
pointing in one direction has the same effect for force anisotropy as the loss of contacts in 
the same direction for fabric anisotropy. As a result, force weakening in the weak network 
is all the more efficient as it leads to lower amount of contact loss. This condition can, for 
example, be achieved for higher level of connectivity, i.e. larger values of z in the critical 
state. 



6 Effect of material parameters 

In this section, we briefly discuss the effect of several material parameters with respect to 
the mechanisms that underly shear strength in granular media. More details will be given 
elsewhere. 

There are several shape parameters that may lead to enhanced shear strength through 
force anisotropy or fabric anisotropy. We consider here polygonal particles as compared 
to circular particles |Azema e t al.(2007)]. The first sample, denoted SI, is composed 
of 14400 regular pentagons of three different diameters: 50% of diameter 2.5 cm, 34% 
of diameter 3.75 cm and 16% of diameter 5 cm. The second sample, denoted S2, is 
composed of 10000 disks with the same polydispersity. The coefficient of friction is 0.4 
between particles and with the walls. At equilibrium, both numerical samples are in 
isotropic stress state. The solid fraction is 0.80 for SI and 0.82 for S2. The isotropic 
samples are subjected to vertical compression by downward displacement of the top wall. 
Figure [7] displays the evolution of a as a function of the cumulative shear strain e q in 
both packings. In both cases, a increases from to its largest value in the critical state. 
Surprisingly, the fabric anisotropy is quite weak in the pentagon packing whereas the 
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Figure 7: Evolution of the anisotropy a with cumulative shear strain e q for a packing of 
pentagons (SI) and a packing of disks (S2). 

disk packing is marked by a much larger anisotropy (~ 0.3). Fig. [8] shows the evolution 
of a n and at- We see that, in contrast to fabric anisotropics, the force anisotropics in the 
pentagon packing are always above those in the disk packing. This means that the aptitude 
of the pentagon packing to develop large force anisotropy and strong force chains is more 
dependent on particle shape than on the buildup of an anisotropic structure. 
According to equation ©, in spite of the weak fabric anisotropy a, the larger force 
anisotropics a„ and a t allow the pentagon packing to achieve higher levels of shear 
strength compared to the disk packing, as shown in Fig. [9] Our numerical data show 
that the strong force anisotropy of the polygon packing results from the edge-to-edge 
contacts that capture most strong force chains, whereas vertex-to-edge contacts belong 
mostly to the weak network. The pentagons provide thus an interesting example where 
the role of fabric anisotropy in shear strength is marginal. Similar conclusions hold for 
polyhedral particles in 3D JAzema et al.(2008)) . 

The effect of the coefficient of friction /.i s between particles on the shear strength was 
discussed in section |2 The saturation of the critical-state friction angle ip c with increas- 
ing /i s is related to the fact that, due to disorder, particle equilibria are fundamentally 
controlled by normal forces. Ideal situations where friction needs to be fully mobilized 
over a large number of contacts exist but are marginal. For example, a column of particles 
each with two contacts may in principal exist, but is of practically zero chance to occur 
within a disordered granular material. The effect of fi s over a c manifests itself through 
Zmin which decreases with /j, s . On the other hand, larger values of /j, s allow for reinforced 
stabilizing effect of weak contacts, increasing thereby force anisotropics and thus shear 
strength. 

The effect of adhesion is to allow for tensile forces mainly in the direction of exten- 
sion between the particles. We find that the tensile forces between particles play the 
same stabilizing role with respect to the strong compressive forces as the weak network 
| Radjai et al.(2001)) . Remark that the privileged direction of weak compressive forces co- 
incides with that of tensile forces. As a result, the main contribution to the shear strength 




Figure 9: Normalized shear stress q/p as a function of cumulative shear strain for the 
samples SI and S2. 



comes from force anisotropy. The fabric anisotropy is generally low and partially inhib- 
ited by the presence of adhesion. Note also that adhesion between particles involves a 
force scale so that its contribution to the shear strength is mainly expressed through the 
Coulomb cohesion c (equation (QJ), but it can also influence the internal angle of friction 
Lp c through fabric anisotropy. 

The size polydispersity is an important factor that affects space-filling properties of gran- 
ular materials. In particular, for a broad size span, the small particles fill and stabilize 
the pores between larger particles. As a result, larger force anisotropies and thus shear 
strengths are expected for higher levels of size polydispersity. Large particles capture 
strong force chains whereas smaller particles are mostly at the center of weak forces 
|Voivret et al.(2008)) . The details of force transmission and force anisotropy depend, 
however, on the size distribution and not only on the span. An expected effect of polydis- 
persity is to allow for higher values of z max and thus enhanced shear strength as predicted 
by equation ( TTOb . 

7 Conclusion 

In this paper, we presented a brief account of physical mechanisms that underly the 
critical-state shear strength of granular materials. The short-comings of the picture of 
granular friction in direct analogy with solid friction was discussed. Recalling the expan- 
sion of the stress tensor in force and fabric anisotropies, a model was presented for the 
accessible geometrical states within a harmonic representation of the microstructure. This 
model, consistent with numerical simulations, relates the critical-state fabric anisotropy 
to two isotropic limit states corresponding to the lowest and highest contact densities of a 
granular packing. The force anisotropy was analyzed in the light of the bimodal character 
of force transmission. It was shown that the shear strength is mainly sustained by the 
strong force network so that force anisotropy is mainly related to the aptitude of a gran- 
ular assembly to build up strong force chains. Finally, the effect of material parameters 
with respect to fabric and force anisotropies was discussed. 
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many useful discussions about granular friction. I present my special thanks to S. Roux 
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